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ABSTRACT 

We carry out a resolution study on the fragmentation boundary of self-gravitating 
discs. We perform three-dimensional Smoothed Particle Hydrodynamics simulations 
of discs to determine whether the critical value of the cooling timescale in units of the 
orbital timescale, /? C rit> converges with increasing resolution. Using particle numbers 
ranging from 31,250 to 16 million (the highest resolution simulations to date) we do 
not find convergence. Instead, fragmentation occurs for longer cooling timescales as 
the resolution is increased. These results suggest that at the very least, the critical 
value of the cooling timescale is longer than previously thought. However, the absence 
of convergence also raises the question of whether or not a critical value exists. In 
light of these results, we caution against using cooling timescale or gravitational stress 
arguments to deduce whether gravitational instability may or may not have been the 
formation mechanism for observed planetary systems. 
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1 INTRODUCTION 

There are two quantities that have historically been used to 
determine whether a self-gravitating disc is likely to frag- 
ment. The first is the stability parameter (Toomre 1964), 



tt£G" 



(1) 



where c s is the sound speed in the disc, Kep IS the epicyclic 
frequency, which for Keplerian discs is approximately equal 
to the angular frequency, Q, E is the surface mass density 



and G is the gravitational constant. Toomre (1964) showed 



that for an innnitesimally thin disc to fragment, the stability 
parameter must be less than a critical value, Qcrit ~ 1. 



|Gammie (2001 ) showed that in addition to the stability 
criterion above, the disc must cool at a fast enough rate. Us- 
ing shearing sheet simulations, he showed that if the cooling 
timescale can be parametrized as 
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u is the specific internal energy and du coo \/dt is the total 
cooling rate, then for fragmentation we require /3 < 3, for a 
ratio of specific heats 7 = 2 (in two dimensions). Otherwise, 
internal heating of the disc due to gravitational instabili- 
ties can stablise the disc. |Rice, Lodato &; Armitage ( 2005 ) 



carried out three-dimensional simulations using a Smoothed 
Particle Hydrodynamics (SPH) code and showed that this 
cooling parameter is dependent on the equation of state. 
They showed that fragmentation can occur if f3 < /5 C rit where 
yftcrit ~ 6 — 7 for discs with 7 = 5/3 and /3 cr it ~ 12 — 13 for 
discs with 7 = 7/5. |Gammie| ( |2001| ) and |Rice et al.| ( |2005 ) 
also showed that in a steady state disc where the domi- 
nant form of heating is that due to gravitational instabili- 
ties, since the gravitational stress in a disc can be linked to 
the cooling timescale in the disc using 
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(4) 



the rapid cooling required for fragmentation can also be 
interpreted as a maximum gravitational stress that a disc 
can support without fragmenting, which they showed to be 
«Gi,max ~ 0.06. On the other hand, other authors have car- 
ried out simulations that suggest that a single critical value 
of aGi,max may not necessarily be the case: |Clarke, Harper-| 
Clark & Lodato (2007) showed that the critical value of ft 
(below which fragmentation will occur if the stability crite- 
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rion is met) may depend on the disc's thermal history: if the 
timescale on which the disc's cooling timescale is decreased 
is slower than the cooling timescale itself (i.e. a gradual de- 
crease in f3) then the critical value may decrease by up to a 
factor of 2. More recently, |Cossins, Lodato &; Clarke (2010) 
showed that the critical value varies with the temperature 
dependence of the cooling law. Given the link between /3 and 
ckgi (equation [4| , this brings aci.max ~ 0.06 into question. 

Further inconsistencies also appear in the literature: 
Rice et al. (2003) found that for a O.1M disc with sur- 
face mass density profile, E oc i? -7 ^ 4 , extending to a ra- 
dius, -Rout = 25 au around a 1M© star, the disc frag- 
ments using /3 = 3 but not for /3 = 5, whereas for a disc 
with mass Mdi sc = O.25M , the disc fragments for f3 = 5. 



2009 ) found that a factor of 2 variation in the particle num- 



On the other hand, Rice et al. (2005) suggested that the 



fragmentation boundary is independent of the disc to star 
mass ratio. iCossins et al. (2010) carried out a simulation 



of a self-gravit ating disc with su rface mass density profile, 



ber did not affect their results. However, these were also for 
non-fragmenting discs. 

In this letter, we carry out a thorough convergence test 
of the value of /3 cr it required for fragmentation using a disc 
set up that is typical of the above studies. Our key result 
is that we are unable to obtain convergence and therefore 
previous results that make conclusions based on a single 
critical value of the cooling timescale (and hence a critical 
value of the gravitational stress) need to be reconsidered. 
We describe the numerical method adopted in Section [2] 
We then describe the simulations and the results in Sec- 
tion [3] Finally, we discuss our results and make conclusions 
in Sections [4] and [5] respectively. 



2 NUMERICAL METHOD 
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(c.f. 



Rice et al. 



(2005 ) who used E oc R 1 ), with Our simulations are carried out using an SPH code originally 



ratio of specific heats, 7 = 5/3, and showed that the criti- 
cal value yftcrit ~ 4. Using equation [4] this is equivalent to 
«Gi,max = 0.1 which also brings the result of aci.max = 0.06 



described above into question. |Meru &; Bate (2010b) recon- 
ciled these inconsistencies and show that the critical value 
of the cooling timescale is dependent on Ei? 2 /M* i.e. the 
ratio of the surface mass density at radius, R, to the stellar 
mass, and concluded that the critical value of the cooling 



time determined by Rice et al. ( 2005 ) is not a general rule 



However, many of their results and the previous results can 
also potentially be explained if the fragmentation is resolu- 
tion dependent. 

Surprisingly, a proper resolution study that demon- 
strates convergence of the fragmentation criteria with in- 
creasing resolution has not been done. Most authors who 
simulate self-gravitating discs consider the resolution crite- 



rion set byjBate fe Burkert| (|1997| (e.g. |Mayer et al.|l2002| 
Lodato fc Ricel|2004| |Stamatellos fe Whitworth||2008| |For-| 



gan et aL]2009| . Some authors carry out a resolution test 



(in addition to ensuring the resolution criterion is satisfied) 
by increasing or decreasing the number of particles by a fac- 
tor of 2 (e.g . |Lodato et aT1|2007| |Clarke et al.|2007| |Cossins 



et al.|[2009| . But large changes in resolution have not been 
tested. In an SPH code (used by many of the above men- 
tioned authors), if the number of particles is increased by a 
factor, /, the smoothing length, and hence the resolution, 
is increased by a factor / d 5 where d is the number of di- 
mensions. |Rice et ah] ([2005 ) carried out a resolution test on 
one of their discs by decreasing the number of particles from 
250,000 to 125,000. They found that the fragmentation re- 
sults appeared to be unaffected by this. However, given that 
the simulations were carried out in three-dimensions, it is 
important to note that this resolution test was equivalent to 
decreasing the linear resolution by only « 21%, so it is un- 
surprising that significant differences were not seen. On the 



other hand, Clarke et al. (2007) carried out a simulation of 



a fragmenting self-gravitating disc and increased the num- 
ber of particles from 250,000 to 500,000 and suggested that 
fragmentation may be affected by resolution. In addition, the 
earlier work by Gammie (2001 ) carried out a resolution test 



for simulations that did not fragment but a resolution test 
was not carried out on the fragmentation boundary. Other 
authors (e.g. Lodato et al.||2007 Cossins, Lodato & Clarke 



developed by Benz (1990) and further developed by Bate, 
Bonnell fe Price| ( |1995| ) "and |Price fe Bate] fl2007| ). We in- 
clude the heating effects in the disc due to work done on the 
gas and artificial viscosity. The cooling in the disc is taken 
into account using the cooling parameter, f3 (equation |2| , 
which cools the gas on a timescale given by equation [3] To 
model the shocks in the discs, we use an artificial viscosity 
(Monaghan & Gingold 1983) with SPH parameters fixed at 
(aspH,/fepH) = (0.1,0.2) (as used by |Rice et al. 2005 and 
Meru fe B ate 2010b). For further details, seelMeru & Batel 
<2010bL 



3 SIMULATIONS & RESULTS 

The disc and star properties used to carry out the simula- 
tions in this letter are exactly the same as those used by 
Rice et al. (2005): a O.1M disc surrounding a 1M© star, 



spanning a radial range, 0.25 < R < 25 au. The initial sur- 
face mass density and temperature profiles are E oc R- 1 and 
T oc i? _1//2 , respectively, and the temperature is normalised 
so that the minimum initial Toomre stability value at the 
outer edge of the disc, Q m in = 2. The discs are modelled 
with a ratio of specific heats, 7 = 5/3. 



Rice et al.|(|2005) modelled this disc setup using 250,000 



particles. In Section 4 of |Meru Sz Bate| (|2010b), we modelled 
this same disc setup with f3 — 5, 5.5, 5.6 and 6 (Simulations 
Benchmark 1-4 of |Meru &; B ate 2010b]) and showed that the 
critical value of the cooling timescale below which fragmen- 
tation occured was /3 cr it ~ 5.6, a value in reasonable agree- 



ment with Rice et al. (2005). We carry out more simulations 



at this resolution, which we summarise in Table [T] 

We also simulate the same disc at higher resolutions by 
using 2 million and 16 million particles (i.e. the smoothing 
length is decreased by a factor of 2 and 4, respectively). In 
addition, lower resolution testing can also be used to test 
for convergence. Therefore, we carry out additional simula- 
tions using 31,250 particles (so that the smoothing length 
is a factor of 2 higher). We simulate the discs using vari- 
ous values of the cooling timescale in units of the orbital 
timescale, /3, to determine the critical value, /3 cr it, at differ- 
ent resolutions. The simulations were run either for at least 6 
outer rotation periods (ORPs) or until the discs fragmented. 
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Simulation name 


No of particles 


P 


Fragmented? 


31k-beta2 


31,250 


2.0 


Yes 


31k-beta2.5 


31,250 


2.5 


Yes 


31k-beta3 


31,250 


3.0 


Yes 


31k-beta3.5 


31,250 


3.5 


No 


31k-beta4 


31,250 


4.0 


No 


250k-beta5 


250,000 


5.0 


Yes 


250k-beta5.5 


250,000 


5.5 


Yes 


250k-beta5.6 


250,000 


5.6 


borderline 


250k-beta6 


250,000 


6.0 


No 


250k-beta6.5 


250,000 


6.5 


No 


250k-beta7 


250,000 


7.0 


No 


250k-beta7.5 


250,000 


7.5 


No 


2m-beta5.5 


2 million 


5.5 


Yes 


2m-beta6 


2 million 


6 


Yes 


2m-beta6.5 


2 million 


6.5 


Yes 


2m-beta7 


2 million 


7 


Yes 


2m-beta8 


2 million 


8 


Yes 


2m-betal0 


2 million 


10 


borderline 


2m-betal0.5 


2 million 


10.5 


borderline 


2m-betall 


2 million 


11 


No 


2m-betal5 


2 million 


15 


No 


16m-betal0 


16 million 


10 


Yes 


16m-betal8 


16 million 


18 


borderline 



Table 1. Table showing the simulations carried out and the key 
fragmenting results. The simulations labelled as borderline are 
those that show fragments forming which then shear apart in less 
than 1 ORP. 



31,250; = 3 250,000; = 5.5 
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Figure 1. Surface mass density rendered images of the fragment- 
ing discs with 31,250, 250,000, 2 million and 16 million particles 
(simulations 31k-beta3, 250k-beta5.5, 2m-beta8 and 16m-betal0, 
respectively). The images are produced at time t = 5.3, 6.4, 5.3 
and 2.5 ORPs, respectively. At higher resolution, the disc can 
fragment for higher values of the cooling timescale. The axes scale 
from -25 au to 25 au in both directions. 



250,000; = 5.6 250,000; = 5.6 




log column density [g/cm 2 ] 



Figure 2. Surface mass density rendered images of the borderline 
cases with 250,000, 2 million and 16 million particles (simulations 
250k-beta5.6, 2m-betal0, 16m-betal8). The left panels shows a 
hint of fragmentation at times, t = 3.8, 4.8 and 6.0 ORPs (top, 
middle and bottom panels, respectively). The right panels shows 
the equivalent simulations a short time later at times, t = 4.2, 
5.8 and 6.2 ORPs (top, middle and bottom panels, respectively). 
Within 1 ORP, the fragments have been sheared apart, classing 
these simulations as borderline. The axes scale from -25 au to 
25 au in both directions. 



Fragments are denned as regions whose surface mass densi- 
ties are at least two orders of magnitude denser than their 
surroundings. Note that there are only two 16 million parti- 
cle calculations as these are extremely time consuming, each 
requiring more than 4 months on 64 compute cores of the 
University of Exeter supercomputer. 

Table [I] summarises the main simulation parameters 
and the key results in terms of whether the discs frag- 
ment or not. The simulations that show no fragments at all 
are classed as non-fragmenting discs. Those discs that frag- 
ment and where the fragments remain without being sheared 
apart are classed as fragmented discs. There are four sim- 
ulations that are classed as borderline (simulations 250k- 
beta5.6, 2m-betal0, 2m-betal0.5 and 16m-betal8). These 
are discs which show signs of fragmentation but the frag- 
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5 6 7 8 

log (number of particles) 

Figure 3. Graph of ft against resolution of the non-fragmenting 
(open squares), fragmenting (solid triangles) and borderline (open 
circles) simulations. The borderline simulations are those that 
fragment but whose fragments are sheared apart and no further 
evidence of fragmentation is seen. The solid black line shows a 
dividing line between the fragmenting and non-fragmenting cases 
and the grey region is where fragmentation can take place. The 
graph shows no evidence of convergence of results with increased 
resolution. The thin dotted line shows how the trend will con- 
tinue if convergence is not reached with higher resolution than 16 
million particles. If convergence can be reached, the dotted line 
is expected to follow a flatter profile. 



ments shear apart rapidly (within 1 ORP) and no further 
signs of subsequent fragmentation are seen. 

Figure [I] shows images of the fragmenting discs in simu- 
lations 31k-beta3, 250k-beta5.5, 2m-beta8 and 16m-betal0 
for f3 values of 3, 5.5, 8 and 10 simulated using 31,250, 
250,000, 2 million and 16 million particles, respectively. Fig- 
ure [2] shows images of the borderline discs in simulations 
250k-beta5.6, 2m-betal0, 16m-betal8. The left panels shows 
the discs appearing to form fragments. However, within 
1 ORP, the fragments shear apart and do not appear to 
form again. 

Figure [3] plots the results in Table [I] The solid black 
line divides the fragmenting and non-fragmenting simula- 
tions and has been included by eye as a fit to where the 
boundary lies. For 250,000, 2 million and 16 million parti- 
cles, the critical value of the cooling timescale in units of the 
orbital timescale, /3 cr it ~ 5.6, 10 and 18, respectively, since 
these have been identified as borderline cases. For the discs 
simulated with 31,250 particles, I assume that the critical 
value is between the fragmenting case of f3 — 3.0 and the 
non- fragmenting case of f3 = 3.5 and thus take /3 cr it ~ 3.25. 
It can clearly be seen that with the data that is available, the 
dividing line between the fragmenting and non fragmenting 
cases increases linearly with linear resolution and therefore a 
convergence has not been reached. Figure [3] also shows how 
the trend may continue if convergence does not take place 
at even higher resolution. 



4 DISCUSSION 

To date, it has been widely accepted that a single critical 
value of the gravitational stress determines if a disc will frag- 
ment. However, a thorough convergence test has not been 
carried out. The dependence of fragmentation on resolution 
presented in this letter shows that convergence has clearly 
not taken place. This therefore calls into question many of 
the previous results on the fragmentation boundary (e.g. 



Gammie|2001j|Rice et al.|2003[ [20051 1 Alexander et al.1 2008 ; 
Cossins et al.||2010| |Meru &; Bate|2010b ), as well as any re- 
sults that have taken the critical value of the gravitational 
stress to be a single value of ckgi ~ 0.06 (or a similarly high 
value of aGi), and based conclusions on these (e.g. Clarke 



2009; Rafikov 2009] |Nero fe Bjorkma n 2009 ; |Kratter et al. 

^olq}. 

It is possible that results that are relative to each other 
still stand. For example, Rice et al. 



may 



2005) clearly 



showed that a stiff equation of state requires a faster cooling 



for fragmentation to occur. Similarly, |Meru fe Bate] ( pOlObfr 
presented physical arguments for why the cooling that is 
required for fragmentation may vary for different surface 
mass densities and star masses. However, until convergence 
is demonstrated, all such results must be treated with cau- 
tion. 

If convergence can be obtained at higher resolution, 
these results imply that /3 cr it can be much higher than pre- 
viously thought. In light of the new results presented here, 
we caution against using cooling timescale or gravitational 
stress arguments to deduce whether certain planetary sys- 
tems may or may not have formed by gravitational insta- 
bility. Clarke (2009) produced an analytical model for the 



structure of a gravitationally unstable disc which is subject 
to realistic cooling. She showed that for optically thick discs 
that are sufficiently low in temperature that they are domi- 
nated by ice grains, 



OiGl 



0.4 



R 



100 au 



(5) 



for a disc with interstellar opacities and surrounding a 1M© 
star, where R is the radius being considered. This relation- 
ship shows that for a maximum value of the gravitational 
stress, a critical radius, i? C rit, can be found outside of which 
fragmentation can occur (for a disc with a shallow surface 
mass density profile). For example, if aGi,max ~ 0.07, then, 
R cr it ~ 68 au. The critical radius also scales as M^ 3 (as 
described by Clarke 2009) and also depends on the metal- 
licity and grain sizes in the disc (as shown by Meru & Bate 
2010al) . One might then use this to compare to observational 



data and conclude whether a system may or may not have 
formed via gravitational instability. 

Table [2] however, shows the equivalent value of QGi.max 
that is associated with the value of /3 cr it identified for each 
resolution considered in this letter, as well as the critical 
radius outside of which fragmentation may take place ac- 



cording to Clarke (2009) for central star masses of 1M© and 
1.5M©. It can be seen that the critical radius of fragmenta- 
tion identified for a disc with 16 million particles is smaller 
than the value identified for a disc with 250,000 by > 20%. 
This can have profound consequences on the interpretation 
of observational results. Table[2]shows that while the value of 
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No of 


Pcvit 


«GI,max 


^crit 


-^crit 


particles 






(Af* = 1M ) 


(M* = 1.5M ) 


31,250 


3.25 


0.12 


77 au 


88 au 


250,000 


5.6 


0.07 


68 au 


78 au 


2 million 


10.0 


0.04 


60 au 


69 au 


16 million 


18.0 


0.02 


53 au 


60 au 



Table 2. Table showing how the critical radius of fragmentation 
according tojclarke (2009) may be affected for a disc surrounding 
a 1 and 1.5Mq star for the different values of /3 C rit identified for 
discs with 31,250, 250,000, 2 million and 16 million particles. 



aci.max from simulations using 250,000 particles may cause 
one to conclude that the outer planet around the 1.5M0 A 
star, HR 8799 ( [Gray fe Kaye|1999[ ), may be difficult to form 
by gravitational instability (since its projected separation 
according to Mar ois et al.| |2008 is « 68 au), the results us- 
ing 16 million particles suggest that this planet could have 
formed by gravitational instability. 

Finally, with the results as they stand, it is possible that 
convergence will never be obtained regardless of the resolu- 
tion. If this is the case, it suggests that the problem may be 
ill posed. In other words, it may not be possible for a disc 
to settle into an equilibrium where there is a balance be- 
tween heating from gravitational instabilities and a simple 
imposed cooling timescale. We notice that once our simula- 
tions have developed significant density structure (but well 
before fragmentation) the temperature of the gas at a given 
radius tends to be lower at higher densities. The cooling rate 
is an imposed function that only depends on radius from the 
star and takes no account of the density structure (which is 
not true of a realistic disc). Therefore, a possible explanation 
for the lack of convergence is that while the low-density gas 
in the disc is being heated by the high-density gas passing 
through it, the high-density gas has no such internal heat- 
ing source and continues to cool forming denser structures. 
The higher the resolution, the worse such a problem may be- 
come. Clearly more testing needs to be done to definitively 
determine the reason for the lack of convergence. 



5 CONCLUSIONS 

We present a resolution test of the critical cooling rate re- 
quired for fragmentation of self- gravitating discs. We find 
no evidence of convergence with increasing resolution. These 
results certainly show that the critical value of the cooling 
timescale is longer than previously thought. However, they 
also open up the possibility that there may be no value of 
/3 for which such a disc can avoid fragmentation, given suf- 
ficient resolution. The latter implies that a self- gravitating 
disc that cools at a rate given by 

du uQ ( . 

dt = T' (6) 

and is only heated by internal dissipation due to gravita- 
tional instabilities, may not be able to attain a self-regulated 
state and will always fragment, regardless of the value of f3. 
This re-opens the question of what the criterion for frag- 
mentation of a self-gravitating disc really is and in addi- 
tion, where in a disc fragmentation can realistically occur. 



These results cast some serious doubts on previous conclu- 
sions concerning fragmentation of self-gravitating discs and 
the results presented here need to be considered when mak- 
ing conclusions as to whether observed planetary systems 
may or may not have formed by gravitational instability. 
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